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This is the third paper in a series describing a numerical implementation of the 
conformal Einstein equation. This paper describes a scheme to calculate (three) 
dimensional data for the conformal field equations from a set of free functions. The 
actual implementation depends on the topology of the spacetime. We discuss the 
\ implementation and exemplary calculations for data leading to spacetimes with one 

spherical null infinity (asymptotically Minkowski) and for data leading to spacetimes 
with two toroidal null infinities (asymptotically A3). We also outline the (technical) 
modifications of the implementation needed to calculate data for spacetimes with two 
. and more spherical null infinities (asymptotically Schwarzschild and asymptotically 

, multiple black holes). 
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In the conformal approach to numerical relativity we give data for the conformal field equa- 
tions on a hyperboloidal initial slice and calculate the conformal spacetime (M, g a b) de- 
termined by the data. The region on which the conformal factor Q, which is one of our 
variables, is positive can be identified with a physical spacetime (M,g a b) satisfying the Ein- 
stein equation. The boundary {Q = 0} of this physical region represents future null infinity 
{J) and possibly future timelike infinity 

Embedding the physical spacetime into a larger conformal spacetime implies a number of 
advantages over conventional approaches: Since null infinity is part of the grid, the deter- 
mination of gravitational radiation is a well-defined and a gauge- ambiguity-free procedure. 
Furthermore we avoid any influence of artificial boundaries onto the result of our calculation, 
and we are enabled to use very efficient high order discretisation techniques, reducing the 
required computational resources by orders of magnitude. 

In the first two papers [I], |2[] in the series we have described the ideas behind the conformal 
approach, the general mathematical background, important properties of the time evolution 
equations, and the code used to integrate the time evolution equations. 
To study properties of asymptotically flat solutions numerically, we need to provide initial 
data, i. e. we have to find solutions of the conformal constraints. In this paper we construct 
data for the conformal field equations by a formally infinite order scheme. The scheme con- 
sists of three basic elements, a Multigrid Newton Method to deal with the non-linearities, 
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pseudo-spectral techniques to achieve formally infinite order, and algebraic multigrid tech- 
niques to invert the linearised elliptic operators. 

In section II we repeat the derivation of the Yamabe equation by the Lichnerowicz ansatz 
and discuss those properties of the solutions which are important for the numerical imple- 
mentation by analytical and numerical means. In the next section we describe the strategy 
for calculating initial data, details of the numerical implementation for the cases leading 
to data for asymptotically Minkowski and asymptotically A3 spacetimes, and give the re- 
sults of two exemplary calculations. In the last section we discuss how the scheme for the 
asymptotically Minkowski data can be extended to data for multiple black hole spacetimes 
by using well-established numerical techniques. 

To avoid too many repetitions, we assume the reader to be familiar with the general ap- 
proach as well as with the equations, both discussed in part I and II of this series ]l], ||, 
and we shall refer to equation (n) of part N by writing (N/n). To the interested reader we 
should also point out the references || [|, ||, |[ [jj which describe a spherical symmetric (ID) 
and an axial symmetric (2D) implementation of the conformal field equations. There is also 
a recent review article by J. Frauendiener in Living Reviews in Relativity ||. 

II. THE YAMABE EQUATION AND THE SMOOTHNESS OF ITS SOLUTIONS 

A permissible set of data for the conformal time evolution equations is given by a set of 
functions which satisfy the conformal constraints. Finding initial data is hence equivalent 
to finding a solution to the conformal constraints. 

The constraints of the conformal field equation (1/14) are regular on the whole conformal 
spacetime (M,g a b), their principal part does not degenerate at J. Therefore, it would be 
nice, if we could solve these constraints directly on our initial slice Et . At present we do 
not know how to do that. The system of conformal constraints is a large coupled system 
of equations and, with the exception of the case of spherical symmetry 0, we do not know 
how to reduce it to a system to which known numerical techniques can be applied, e. g. a 
system of elliptic equations. 

The way we construct initial data for the conformal field equations is an indirect one, we first 
solve the constraints of the vacuum Einstein equations and then construct the rest of the 
data from the conformal constraints. This amounts to providing a numerical implementation 
of the approach which has been used in |9|] by L. Andersson, P. Chrusciel, and H. Friedrich 
to prove existence of regular data for the conformal field equations. With this ansatz we 
are led to solving an elliptic equation of second order. A generalisation of the ansatz in || 
leads to a system of four coupled equations which has been analysed in |[21|| . We will restrict 
ourselves to the case of the ansatz of || for simplicity. In this case we already have the 
freedom to prescribe the conformal metric on the initial slice. Although this is not the full 
parameter space for the initial data, we expect to have covered enough of the parameter 
space to study interesting phenomenae. 

A. The Yamabe equation 

Let M be an asymptotically flat spacetime and £ a hyperboloidal submanifold, i. e. a 
spacelike hypersurface extending to future null infinity. We denote by h a b the induced 3- 
metric, by k a \, the induced second fundamental form, and by k its trace. The Hamiltonian 
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constraint then reads 

(3) R - k ab k ab + k 2 = 0, (1) 

where {3) R is the Ricci scalar associated with h ab . If we denote by (3) V a the covariant 
derivative induced by h ab , the vector constraint reads 

{3 V b k ab - i3 V a k = 0. (2) 

We want to construct initial data representing this geometric situation. To simplify our 
calculation we assume that on our slice 

k a b = l^habk- (3a) 



The vector constraint (0) implies 
We assume 



k = const ^ 0. (3b) 



k > 0. (4) 



Think of S as being smoothly embedded into the conformal extension of our spacetime 
through future null infinity J + and denote by S the closure of E in this extension and by 
S the boundary of S. Let Q be a boundary defining function with non-vanishing gradient 



on S, positive in the interior |22| and vanishing on S, and let h ab be an (almost) arbritrary 
metric on S. 

To reduce the Hamiltonian constraint (|l|) to an elliptic equation for a scalar function 0, we 
make the so-called Lichnerowicz ansatz, which reads in our case 
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1S\ h ab . (5) 



With this ansatz h ab is singular at S indicating that S represents an infinity. 
The Hamiltonian constraint becomes the so-called Yamabe equation, 

4 Q 2 (3) A0 - 4 n( (3 V a 0) ( (3 V a 0) 

- Q (3) R n 2 + 2n (3) Afi - 3( (3 V a fi) ( (3 V a fi)) - h 2 (j) 5 = 0, (6) 

where (3 V a , (3) A, and {3) R are the covariant derivative, the Laplace operator, and the Ricci 
scalar associated with h ab . The Yamabe equation is the equation, which we are going to 
solve numerically. It is an "elliptic" equation with a principal part which vanishes at the 
boundary. 

After having solved the Yamabe equation (§) the remaining members of a minimal set of 
data which consists of (hab, k a b, ^o) (cf. II) are given by 

kab = \kh ab (8) 
n = i (pk - k) , (9) 
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where k is an arbritrary function. Its choice is pure conformal gauge, it does not influence 



the physical spacetime obtained (cf. also [10, section 3] or [11, subsection 2.1.1]). 

To calculate a complete set of data for the conformal field equations from a minimal set of 

data we use the conformal constraints to calculate 7 a h c , (0,1) R a , CM) -Ra a > ^a> <^> f^^ftab : = 

Q <1A) R ab , f Ba b '■= ^B ab , and f Eab = QE ab - 1/2 (1 ' 1] R ab (cf. II for more details). 

Due to our assumption (|3a]) the magnetic part of the Weyl tensor B ab vanishes. For the case 

of vanishing conformal extrinsic curvature, i. e. k — 0, the quantity (0,1) -R a also vanishes. 

The following theorem by L. Andersson, P. Chrusciel and H. Friedrich || gives existence 

and uniqueness of positive solutions to the Yamabe equation and guarantees regularity of 

^R ab and E ab : 

Theorem 1 (Andersson, Chrusciel, Friedrich) Suppose (E to , h ab ,Cl) is a 3- 
dimensional, orientable, compact, smooth Riemann space with boundary S, Q = on 
S and positive elsewhere. Then there exists a unique positive solution <fi of equation (^) and 
the following conditions are equivalent: 

1. The function cf) as well as the corresponding complete set of data, determined on the 
interior £ io ofE to , extend smoothly to all of (S to , h ab ) . 

2. The electric part of the conformal Weyl tensor QE ab goes to zero at S. 

3. The extrinsic 2-curvature induced by h ab on S is pure trace. 

If we did not require, that the extrinsic 2-curvature of S were pure trace, f(i,i)fi ab and fEab 
would not be vanishing on the boundary and therefore {1,r) R ab and E ab would not be regular 
there. For this reason we choose the conformal 3-metric in such a way that the tracefree 
part of the extrinsic 2-curvature of S with respect to h ab vanishes. 

Observe that if we want to have a regular positive solution <p of (|B|) we cannot specify 
boundary values for <p on S. Since for a regular solution the principle part vanishes at the 
boundary, we must have 

(V a fi)(V a fi)-^V 
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0. (10) 



B. Smoothness of solutions of the Yamabe equation on the initial slice 

In our setup S to will be a true subset of the initial slice E 4o represented by the grid. To 
construct a complete set of data on E <0 one is tempted to give fl and h ab on T, to with Q < 
outside S io and solve (|G]) on E to . In the general case we obtain a complete set of data which 
is not continuous on S. Since our discretisation of the time integrator requires sufficiently 
smooth data we cannot proceed this way. In the following we will give an argument, why 
we have to expect a non-continuity, and we will show a numerical example. 

To calculate the curvature variables {1,1) R ab and E ab from /(i,i)i? a 6 anc ^ fEab we have to divide 
by Q. To calculate the limit for Q — > 0, which is the value at S, we apply rHopital's rule 
and get 

/ 

n 



n=0 ^V a n ( 3 >V a ft 
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For the highest derivatives of which appear in the expressions of (1,1) i? a b and E ab we find 

^R ab ~ <9 3 0, (12a) 
E ab ~d 4 (P. (12b) 

Let us now look at the following simple situation of figure |T[ Suppose we have given a 




FIG. 1: A Kruskal-Schwarzschild slice en- 
closing a Minkowski slice. 



spherical symmetric 0, a spherical symmetric h ab such that <Sj is the unit sphere, and k = 0. 
Suppose is positive in the part labeled "Minkowski", which has one boundary diffeo- 
morphic to S 2 , and negative in the part labeled "Schwarzschild", which has two boundary 
components of topology S 2 . Then, due to Birkhoff's theorem, {Cl > 0} is a hyperboloidal 
slice of Minkowski spacetime, and {fl < 0} is a hyperboloidal slice of a Schwarzschild space- 
time, since there are no spacelike slices in Minkowski space connecting two null infinities. 
The latter must have mass m^fl, whereas the mass of the Minkowski part necessarily van- 
ishes. The Bondi mass of the initial slice is given by the integral over S with an integrand 



which is a polynomial expression of our variables involving components of E ab | 12| . Since 
the mass is 0, if we take the integral on the interior side, and non-vanishing, if we take the 
integral on the exterior side, the integrand cannot be continuous. The only source of the 
non-continuity can be the derivatives of 0. From the proof of theorem p] it follows that is 
at least C 3 across the boundary S. Since the highest derivatives of in E ab are of fourth 
order, can in general not be C 4 across the boundary S. 

In a numerical code we would, of course, give boundary values for on the grid boundaries 
instead of specifying where the outer S is placed. However, unless we happen to prescribe 
by coincidence the appropriate boundary data on the grid boundary, we can expect to get 
the same behaviour. To check this hypothesis we have performed numerical test calculations 
for asymptotically A3 data with one Killing vector (a 2D calculation). 
In the following we give the results of a typical numerical calculation. In this 2D calculation 
we have given the same free functions as in [J], 

ti = Ul-x 2 ) (13a) 
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h xx = 2e- 2xcos{y) 
h xy = 2ti (x 2 — sin(?/ 2 )) e" 

h yy = 2 (l + Q 2 (x 2 - sm(y 2 

h zz = 2 e - 2xcos{y) 
k = 0, 



x cos(y) 



(13b) 
(13c) 

(13d) 

(13e) 
(13f) 



where (x, y) G [—1.4, 1.4] x [— tt, tt]. The free functions satisfy the regularity condition on S. 
We then solve the Yamabe equation by discretising it by centered second order stencils, 
inverting the resulting sparse matrix with the AMG library |]13, 14]], and taking care of the 



non-linearities by a Multigrid Newton Method [ |i~5[ 
Although the free functions depend on both coordinates, we plot the results only for y = 
— 7r/4 for the reason of simpler graphics. The solution eft (solid line) for the calculation with 
the highest resolution performed, which is a 640 x 640 grid calculation, is shown in figure |2|. 
We have also plotted the conformal factor Q (dashed line) to indicate the location of the Ss 
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FIG. 2: Plot of (ft (solid line) and (dashed line) along the 
y = 7r/4 line. 



at x = — 1 and x = 1. The function (ft seems to be perfectly smooth across the Ss. 

By calculating the quantities {1,1) R a b and E a b we see that this is not the case. To do so we use 

the second order scheme described in II. Figure |3| shows for the grid sizes of 20 x 20, 

40 x 40, 80 x 80, 160 x 160, 320 x 320, and 640 x 640 gridpoints (by (M) i?^ we denote the 
x coordinate components of the tensor (1,1) i? a b). Increasing line density correlates with finer 
grids. Obviously convergence near the two Ss is slow — the difference between the various 
grids is larger — and obviously {1,1) Rxx_ converges against a function which is only C° at the 
two Ss. This non-smoothness also explains the slow convergence rate. 
Convergence is even slower in figure [|, where we plot E m . The solution slowly converges to 
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FIG. 3: {1A) Rxx for grids from 20 x 20 (dotted line) to 640 x 
640 (solid line). 



a function with steps at the two Ss as expected from the theoretical considerations in the 
first part of the subsection. 

III. CALCULATING EXTENDED HYPERBOLOIDAL INITIAL DATA THE 

NUMERICAL IMPLEMENTATION 

A. A rough description of the main ideas 

What has been said in the previous section suggests the following strategy to calculate 
extended hyperboloidal initial data: 

1. Solve the Yamabe equation (H) on S to to get a minimal set of data. 

2. Calculate a complete set of data on E to . 

3. Smoothly extend the data from S io to S io . 

The main part of the numerical implementation of step [l] and step |2| is solving "elliptic" 
equations with a principal part which degenerates at the boundary. To be able to obtain 
sufficient accuracy with an acceptable consumption of computer resources even without 
symmetry assumptions (3D) we use pseudo-spectral methods similar to the 2D calculations 
described in |7j . Since there exists a huge amount of literature on pseudo-spectral methods, 
we only very briefly describe our choices. The reader who wants to learn about the basics is 



referred to the literature, e. g. |T6[ for the general theory or [17j for other general relativistic 
applications. 

Since pseudo-spectral methods are formally infinite order, the grid on which we discretise 
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FIG. 4: JS^ for grids from 20 x 20 (dotted line) to 640 x 640 
(solid line). 



and invert the elliptic operator can be pretty coarse. The accuracy achieved is very high, 
as soon as the spectral grid resolution is sufficient to represent the structure of the solution. 
The limit of achievable accuracy is given by the accumulation of rounding errors in the fast 
Fourier transformation, the matrix inversion, and the consecutive solving of elliptic equa- 
tions. Typically, using 32 gridpoints in each dimension is sufficient to achieve an accuracy 
which cannot be further improved. 

For step ||] we use the representation of our variables in the spectral basis to calculate a 
smooth extension. 



B. The Yamabe solver and the f/Cl divider 
For simplicity we write the Yamabe equation in the form 

E 0-lp0 5 = O, (14) 

where E is a linear operator. 

To deal with the nonlinearity we use the Multigrid Newton Method for nonlinear systems 
described in [p~5| 



Assuming we have given an approximate solution 0^ we calculate a hopefully better solution 

0(n+l) = 0(n) + e ( 15 ) 

where the correction is the solution of 

EM := E5<f) - 5 (> (n) ) 4 5<j) = -E0 (n) + (0 (n) ) 5 , (16) 
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which is obtained by linearising the nonlinear equation (0). Since 5(f) is a correction to the 
solution, the order of the scheme and the eventual quality of the approximation depend on 

the order and the accuracy of the calculation of the residuum — E0( n ) + ^0*™^ , but not on 
the discretisation of E/. 

The value of the parameter e is determined in the Multigrid Newton Method 
following algorithm: We take the first value of e in the sequence 1, |, |, . . . for which 




for an appropriate norm ||.||, in our case the L 2 norm. This adaption of e on each step 
stabilises the scheme and makes it less sensitive to a good initial guess Numerical 
experiments with a fixed e were only stable for e < l/2 n , where n is the dimension of the 
hypersurface S to . 

The residuum is calculated by pseudo-spectral methods, which means, that from the values 
of (p^ at the gridpoints we calculate the spectral coefficients. Then we differentiate in 
the spectral space, where differentiation is essentially a multiplication, and transform back 
to the grid to get high order approximations of the derivatives. With these values of the 
derivatives we calculate the value of the residuum on the gridpoints. 

The order of the differentiation with respect to a coordinate x a is approximately the number 
of functions in the spectral basis for this coordinate, which is approximately the number 
of gridpoints in the corresponding dimension (a more accurate statement depends on the 
spectral basis, see below). Therefore the order of the differentiation increases with the 
number of gridpoints. This increase of the order with the number of gridpoints is called 
"formally infinite order". 

For reasons of numerical efficiency the choice of the spectral basis is restricted to those 
bases for which we can use fast Fourier transformation (FFT) techniques to perform the 
transformation to and from the spectral space. 

To perform the FFTs we use the FFTW library by M. Frigo and S. G. Johnson, which is 
a C subroutine library for computing the discrete Fourier transform (DFT) in one or more 
dimensions, of both real and complex data, and of arbitrary input size [jnj. The library 
includes a parallel version for POSIX threads and the FFT is done by an NlogN algorithm 
even if the prime factor decomposition of N contains large prime numbers. This has large 
advantages over simple NlogN algorithms which are based on the assumption that the 
number of gridpoints N is a power of 2, e. g. when it comes to the issue of grid refinements 
in 3D. 

The derivatives in the operator Ei are discretised by three point stencils. These stencils 
depend on the spacing of the grid which is given by the choice of the spectral basis. We are 
going to describe it below when we discuss the specific cases. Due to the dependence on the 
spectral basis we call this grid a spectral grid. In general the spectral grid does not coincide 
with the grid used for time evolution and is significantly coarser. 

Where the boundary of the spectral grid coincides with S, the initial approximation (f>^ for 
the solution is chosen such that it coincides with on S. Then the boundary value for 
5(f) is in each iteration. If the grid boundary is inside physical spacetime, which is the case 
for more complicated forms of E to , where we need overlapping spectral grids to cover the 
physical part of the initial slice (multiple black holes, subsection |TV|), the boundary values 
for 5(f) are deduced from the other grids and iterated in a Schwarz alternating procedure. 
After having provided boundary values we use algebraic multigrid techniques to solve for 5(f). 
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To do so we use the AMG library ]13|, [14|], which was kindly provided to us by K. Stiiben from 
the Gesellschaft fur Mathematik und Datenverarbeitung. Different topologies of the spaces 
on which we solve the elliptic equations show up in different structures of the discretisation 
matrix of Ej. Since AMG automatically derives a multigrid coarsening strategy from the 
structure of the discretisation matrix, only minor changes need to be made in the code which 
inverts the elliptic operator to adapt it to the various cases described below. This is a huge 
advantage. The price to pay is the fact, that the AMG library is the only part of the code 
which does not give excellent parallel scaling, in fact the AMG version used does not run in 
parallel at all. For a typical 3D run, which produces data for a time evolution grid of 100 3 
or more gridpoints, the total time spent in the AMG library is of order of a few percent of 
the time spent in the rest of the code, although the rest has excellent parallel performance. 
Therefore, the non-scaling of the AMG part of the code is not a serious issue. 



To calculate g = f /Q we use the same numerical techniques applied to equation (II/7), 
which can formally be written as 

Gg = F[f] (18) 

and has the "linearised" form 

G l 6g = -GgW+F[f\. (19) 

Observe that the metric used in (II / 7) does not need to be identical with the 3-metric h^. 
And indeed, the code becomes much simpler if we use the diagonal Euclidean 3-metric 5 a b- 
As boundary values we use (II /8) on the grid boundaries which coincide with S and f/fl 
on the grid boundaries which lie in the physical part. Of course, due to the latter, the grid 
boundaries which lie in the physical part must not intersect S. As initial approximation we 
use a low order approximation of equation (II / 7) . 

The attentive reader may ask the following. Equation (II/7) is a linear equation for g, what 
is the reason for doing an iteration which is typically used to deal with nonlinear terms. 
The reason is very simple. If we discretised G to low order, we would get a sparse discreti- 
sation matrix, but also an inaccurate solution. If we discretised with high order stencils, 
we would get a not- very-sparse discretisation matrix, which is difficult to invert, and a 
complicated boundary treatment. By using the iteration on the linear equation we get a 
formally infinite order solution, which turns out to be very accurate, but we still have a 
sparse discretisation matrix. The repeated inversion of G/ is the price to pay. 



C. The case of asymptotically A3 spacetimes 

From the viewpoint of the initial data solver the_ simplest class of initial data are asymptot- 
ically A3 spacetimes. There the physical part S to of an initial slice has topology R x T 2 , 
and the boundary S consists of two two-dimensional tori T 2 (figure 4 of I). By doing an 
appropriate coordinate transformation within the initial data surface we can always ensure, 



that the coordinates y and z are 2n periodic and parametrise the two-dimensional torus p3 
and that the two parts of the boundary S, in the following called Si and S r , are at constant 
x values x = a and x = b. Although the code allows arbritrary values of a and b, as well as 
other periodicities in the y and z direction, we only describe the case with a = — 1, b = 1, 
and 2ir periodicity for simplicity. 

Since y and z are 27r periodic we can simply use the Fourier functions e tkyV and e lkz 2 as basis 
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in the y and z direction. The spacing of the spectral grid in y and z direction is equidistant, 
the Ny jZ nodes of the spectral grid are placed at 

VJ ■= ^ (20a) 

2-rck . , . 

:= -TT- ( 20b ) 



N. 



z 



In the x direction we use Chebychev polynomials 

Tkx( x ) '■= cos (kx arccos(x)) . (21) 
The gridpoints are placed at the Gauss-Lobatto nodes 

Xi -= cos (22) 

where N x + 1 is the number of gridpoints in the x direction. 
For every function 

/:[-M]-i2, (23) 
there is a one-to-one relation to the function 

/ o cos : [0, tt] -> i2. (24) 

We can reduce the Chebychev transformation to a Fourier transformation, which of course 
can be treated by FFT techniques, by transferring the non-equidistant gridpoints in [—1, 1] 
to equidistant gridpoints in [0, tt] and extending the range to the interval [0, 2 r w[ by applying 
the symmetry 

cos(7r + x) = — cos(x). (25) 
For a function f(x, y, z) the spectral representation 

N x N y /2 N z /2 

f(x,y,z):=Yl E E °>k xkykz T kx (x)e ikvy e ik * Z (26) 

k X =0 ky = ~Ny/2 k z = -N z /2 

is given by the coefficients a kxkykz . Depending on N y and N z some of the coefficients are 
zero and due to the reality of / the complex numbers a kx k y k z are not independent. 
Then the coefficients a xyz kxkykz of the derivatives d xyz are given by 

a v k xky k z = k y a kxkykz (27a) 
a z k x k y k z = k z a kxkykz (27b) 

and the recursion relation 

| d X N X kyk Z = ^N^l kyk Z = ^ X 

\ck x a x k x k y k z = a x k x +2k y k z +'2{kx + ^)a x k x +ik y k z , k x = N x — 1, . . . , 0, 
where Co = 2 and c kx — 1 for k x > 1. 

Due to ( |27cj ) we loose one discretisation order with each d x . 
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There exist FFT algorithms which take into account the even symmetries of the func- 



tion (23), the so-called fast cosine transformations [19]. They are significantly more com- 
plicated than a normal FFT, but reduce the required space by a factor of 2 approximately. 
Since the space temporarily required in the initial data solver part of the code is significantly 
less than the space temporarily required in the time evolution part, we abstained from mak- 
ing use of the fast cosine transformation techniques and the implied memory savings. 
Therefore, to calculate the spectral representation on the spectral grid, we proceed as 
follows: We first relabel the grid covering [—1,1] x [0, 27r[x [0, 2ir[ to equidistantly cover 
[0, 7r] x [0, 27r[x [0, 2ir[ and then extend it to [0, 2tc[x [0, 2tc[x [0, 2tc[. We use a three dimen- 
sional FFT for real numbers on the extended grid to calculate the spectral representation. Al- 
though our scheme can deal with any number of spectral gridpoints (N x , N y , N z ) > (2, 2, 2), 
for reasons of numerical efficiency it is advisable to use values for which we get efficient 
FFTs. Of course, the large grid covering [0, 2n[x [0, 2n[x [0, 2tt[ is only used for calculating 
the spectral transformation and its inverse, tasks like calculating the residuum and the in- 
version of the discretisation matrix of E; are performed on the original, smaller grid covering 
[—1, 1] x [0, 27r[x [0, 2ir[. The procedure for the inverse transformation is obvious from what 
has been said. 

To derive the discretisation matrix of E; we discretise the derivatives in E; as follows 

°W — > ~~, r~7 ~Ji-i,j,k H 7 r~, r~Ji,j,k 

— %i-l) \ x i — %i-l) ~ %i) \%i ~ 

+ 7 " T7~^ sfi+l,j,k ( 28a ) 

— %i-\) — %i) 

9 y f —> fi,j-i,k H fi,j+i,k (28b) 

Vj+i - Vj-i Vo ■ i - Vj-i 

&zf ~^ fi,j,k-i H fi,j,k+i (28c) 

Zk+l — z k-l z k+l — z k-l 

<^xf ~~ ^ 7 w ~fi-l,j,k ~ ~ ( w \fi,j,k 

{%i+l — {%i — — Xi) \%i ~ 

+ 7 \ ( (28d) 

[Xi+l — Xi-l) [ x i+l — %i) 
9yf —> 7 72 fij-hk - 7 72/iJ-fc + 7 72 fi,j+hk ( 28e ) 

[Vj+i - Vj-i) [Vj+i - Vj-i) [Vj+i - Vj-i) 

2 4 8 4 

®zf ~^ 7 72fi,j,k-l - Z^fiJ.k + 7 72 fi,j,k+l ■ ( 28f ) 

[Zk+l — Zk-l) {Zk+l — Zk-l) [Zk+l — Zk-l) 

The mixed second derivatives are the obvious combinations of the first derivatives. 
Combining the previous with subsection |l 1B| we have a very efficient elliptic solver for the 



Yamabe equation (cf. subsection [III E] for numbers). The solution of the Yamabe equation 
determines a minimal set of data. 

Next we calculate the first and second derivatives of the minimal set of data by pseudo- 
spectral methods, which immediately leads to values for +- a b c , ((U) i?a, {l ' 1) R a a , &a, 
f(i,i)R a b := ^ <1,1) ^afe, and fEab = QE ab — 1/2 (1,1) Rab on the spectral grid. Due to assump- 
tion (|3a|) fsab '■= QB a b is identical 0. 

To calculate CM) i£ 6 and E a b on the spectral grid we apply the procedure for dividing by Q 
as outlined in the second part of subsection |III B| . 

We now have a complete set of data on the spectral grid for S io . We need to transform these 
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data to the grid which is used in the time evolution and which extends beyond S io . Analyt- 
ically we could just evaluate the sum in ( p6|) for any (x, y, z) on the grid. Numerically that 
is not possible: For |x| > 1 the absolute value of the nth Chebychev polynomial grows as 
fast as \x\ n . Although the Chebychev coefficients of analytic functions decay exponentially 
fast for sufficiently large n, the numerically calculated coefficients do not decay to exactly 
due to rounding errors, typically not larger than 10~ n , but non- vanishing. Due to the rapid 
growth of \x\ n for large n, these rounding errors would dominate the numerical result. The 
described way of extending is therefore numerically unstable. 

A possible solution goes as follows: For \x\ > 1 we replace the Chebychev polynomials Tk(x) 
by some functions Tk(x) which fuse sufficiently smooth into the Chebychev polynomials at 
| x | = 1 and which have bounded growth. There are of course an unlimited number of options 
to do so, or choice is 

rf, / x f (sigm;) fc cosh ( k cosh _1 x) \x\ > 1 

Tk{x) = \ K ° J \ J (29) 

[ cos {k cos x) \x\ < 1 

with x = (tan ^ (x — signx)) / (jj-j + signx. We call the process of calculating initial values 
on the time evolution grid from the spectral grid "extension procedure" . 



D. The case of spacetimes which are asymptotically Minkowski 

Asymptotically Minkowski spacetimes are spacetimes whose null infinity has spherical cuts. 
Therefore the topology of St is the one of a three-dimensional ball (cf. figure 1 of I). 
We assume that (x a ) = (x,y,z) is a coordinate system induced by R 3 , which we call a 
Cartesian coordinate system, and that the Cartesian components hgb of the 3-metric h a b 
are smooth. These are the coordinates in which the time evolution is done. We call the 
coordinate system (r, ip) defined by 

x = r sin $ cos ip (30a) 
y = r sin t? sin ip (30b) 
z = r cost? (30c) 

the polar coordinate system x a> . Without loss of generality we can assume that S coincides 
with the r = 1 coordinate surface. The interval (r, ip) e [0, 1] x [0, n] x [0, 2tt[ then covers 
the unit ball B 3 . 

Using polar coordinates has the advantage, that the boundary S coincides with a coordinate 
surface. On the other side we have to deal with the coordinate singularities at d = or tc 
(the z-axis) and r = (the origin). The coordinate singularities cause singular behaviour of 
the metric and its derivatives and potentially lead to numerical instabilities. To avoid these, 
a careful choice of the placement of the grid nodes and special measures in the spectral 
transformation must be taken. In this paper we describe what we have done. The reader 
can find a more exhaustive discussion of what can be done in . 



Apart from the changes enforced by the coordinate singularities, the code is very similar to 
the previous case of asymptotically A3 data. 



Since any interval [a, b] can easily be mapped to [—1,1], functions on any interval [a, b] can 
be represented by Chebychev polynomials. Therefore, one could in principle work on the 
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r-interval [0, 1]. Due to the coordinate singularities, this would not be a good idea. Firstly, 
there would be a Gauss-Lobatto node directly at the origin. And secondly, near the origin the 
gridpoints would be extremely dense. Instead of representing the ball by [0, 1] x [0, it] x [0, 2tt[ 
we use [—1, 1] x [0, tt[x [0, 7r[ as minimal coordinate patch. By the symmetries 

f(r,# + Tr,<p) = ±f(-r,#,<p) (31a) 
f(r, 0, <p + tt) = ±/(-r, n-#,<p), (31b) 

where the signs depend on what tensor component is represented by /, the use of the 
function (p4"|) instead of (0), and the identity (p5|), this interval can be mapped to the 
interval [0, 27r[x [0, 27r[x [0, 2ir[, on which everything is 2ir periodic. Spectral transformations 
are performed on this extended range [0, 27r[x [0, 27r[x [0, 2ir[. As in the previous section the 
use of the extended range is a waste of memory, but again this is not influencing the total 
memory requirements of the code. 

Now we have to place the gridpoints in such a way that there is no gridpoint at the origin 
and on the axis. We use N r + 1 gridpoints, where N r is odd, to cover the r-interval [—1, 1]: 

Ti = cos(z— ) i = 0, . . . , N r . (32a) 

Since N r is odd, there is no gridpoint at r = 0. The interval [0, 2tt[ of the extended range is 
then covered by 2N r gridpoints, i. e. 4 is not a divider of 2N r , and FFTs based on powers 
of 2 cannot be used. 

To cover the ^-interval [0, ir[ we use N# gridpoints which we place at 

^=j-^ + -^-, j=0,...,iV,-l. (32b) 

The shift of ensures, that there is no gridpoint at the axis. In total we have 2N$ 
gridpoints on the ^-interval [0, 2tt[. 

To cover the (^-interval [0, tt[ we use N v gridpoints which we place at 

<Pk = k^-, k = 0,...,Ny-l. (32c) 

There is no shift needed, any value of if can be assumed. In total we have 2N ip gridpoints 
on the yj-interval [0, 27r[. 

The discretisation of the derivative operators in Ei happens as in (28), where (x,y,z) is 
replaced by (r, </?). 

We now have a grid which avoids the coordinate singularities. This is still not sufficient 
to avoid numerical problems. The problems are caused by the components which behave 
like powers of 1/ sin($) near the axis and/or like powers of 1/r near the origin. If we make 
a spectral transformation we always get significant values for the highest frequencies, no 
matter how many basis functions (=gridpoints) we use. These high frequencies make the 
Multigrid Newton Method unstable. To avoid the high frequencies we only apply spectral 
transformation to "regularised quantities". What we mean by "regularised quantities" will 
become obvious from what follows. 

As free functions we give the Cartesian components hgb of the 3-metric h ab , Cl, and k as 
functions of the polar coordinates. We suppose that the Cartesian components of the metric 
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and k are given in a form which extends to the whole time evolution grid. The polar 
coordinate components h a i b > of the 3-metric h a b are then given by 

h rr = (2c {p c^h xz + 2c$h yz s ip ) s$ + h zz (c^) 2 + [^c^hxyS^ 

hrd = r Cphx^cv) 2 + h yz s^(ci}) 2 + s$ (—(c#h zz ) + 2c ip c#h xy s ip + c^h^c^) 2 + c^h yy (s ip y 



rh 



h 



rtp 



'Ctphxx^ip CiphyyS^ -\- hxyictp) hxyi^s^ J 



hm = r 2 (-2cpCoh xz - 2c^h yz s (p ) s# + 2c (f h xy s (f (c^) 2 + h X x(c<p) 2 (c & ) 2 

+ hyy(Ci)) 2 {S V ) 2 + h ZZ (S, & ) 2] 

1&& 



r 2 L 



(33b) 
(33c) 

(33d) 



J Si) 



r s$h 



2 2 

r s# 



dtp 

2c l ph X yS ! p ~t~ hyyi^c^ ~\- h X x(s t 



°xx \°tp J 



r 2 (so) 2 h 



(33e) 



(33f) 



where c# := cos$, c v := cosy?, s$ := sin??, s v := sin<£>, and we have written h rr although we 
really mean h rr . The inverse metric h ab can then be written as 



a'b' 



( h rr 



•sin'!? _ 



[ -^—M Tl c t4-5^ T-^sh 



r 2 sini? 

1 z,<w 



(34) 



where h-- is the inverse of the matrix h a >y. 

To calculate the polar coordinate components 7-&V °f the 3-Christoffel symbols 7 a & c and the 
3-Ricci scalar (3) i?, which contains derivatives of the 7- b > c rS, we first calculate the first and 
second derivatives of the "regularised quantities" h a i b > and h-- by spectral techniques, and 
then put in the singular terms by hand. This avoids the high frequency problem mentioned 
above. 

The solution of the Yamabe equation (|6|) is then calculated as a function of the polar 
coordinates as before. Since the scalar is regular on the axis and at the center, no regu- 
larisation is needed. 



From the solution we calculate a complete set of conformal data as follows. The Cartesian 
components of h ab are taken from the free functions evaluated at the gridpoints. Then, the 
Cartesian components of the 7 a b c s are obtained by fourth order differentiation (II /10a) on 
the time evolution grid. We do not use the extension procedure for the 7 a b c s for two reasons. 
Firstly, the extension procedure is slow. And secondly, transforming the polar coordinate 
components of the 7 a b c s to Cartesian components requires a careful treatment of the regu- 
larity conditions. As price for this convenience we get only 4th order approximations of the 
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Cartesian components of 7 a 6 C on the time evolution grid instead of infinite order approxi- 
mations. But even if we started with infinite order data, after some time steps everything 
would only be 4th order, since the time evolution scheme is "only" 4th order. 
Then we use the extrapolation procedure to calculate Q on the whole grid. The spatial 
derivatives Q a of Q are again obtained by fourth order differentiation. The time derivative 
Qo of Q is calculated from the free function on the time evolution grid directly. 
The scalar ui is calculated on the spectral grid and then extended. The quantities (0,1) R a , 
/(x,i)r , and fE ab are tensor components. We first calculate the polar components on the 
spectral grid and then transform to Cartesian components. To the Cartesian components 
of (01) -R a we apply the extension procedure. 

The Q division is applied to the Cartesian components of f(i,i)R an d fE ab on the spec- 
tral grid. This leads to the Cartesian components of (1 ' 1] Rab and E^. We complete the 
construction of the complete set of initial data by extending those to the time evolution 
grid. 



E. Two examples 

In this subsection we give two examples for initial data calculated with the solvers described 
above. The first example is an asymptotically A3 data set, the second an asymptotically 
Minkowski data set. We have given very general free functions. The purpose is to demon- 
strate the performance of the code described. In a realistic parameter study one would 
probably choose much simpler free functions. 

For the asymptotically A3 case we have chosen 

Q = ~ (l - x 2 ) (35a) 

/ l + ti 2 (smy) 2 ifi 2 (x 2 -(siny) 2 ) \Cl 2 (cosz) 2 \ 

hab= \Q 2 {x 2 - {siny) 2 ) l + fi 2 (sin^) 2 |fi 2 cos(x7r) (35b) 

V ±fi 2 (cosz) 2 |fi 2 cos(x7r) 1 + |^ 2 / 

k = cosy. (35c) 

Without the fr terms, the choice of leads to data for an A3 spacetime. 

Figure [5| shows the convergence of the violation of the constraint (3) V&-Ex b + {3) exbck bd Bd c = 

(II/14d, with a = x), which is the constraint which shows the largest violation, for increasing 

spectral grid density. For simplicity we have plotted the L 2 norm over (y, z) as defined in 

(11/30). We observe that the violation is rapidly going to as could be expected. 

As we have 90 constraints and the violation of (II/14d,a = x) dominates for all spectral 

resolutions, the average violation of constraints is more than by a factor of 50 less. 

For these runs we stopped the Multigrid Newton Method iteration when the L 2 residuum 

dropped below 10~ 7 , at that time the approximation changed only by an order of 10 -11 per 

iteration. To achieve such a small residuum we typically need less than 20 Iterations, which 

means that the average improvement of the residuum per iteration is a factor of 3. 

For the 33 x 32 x 32 spectral grid the average computer time per AMG step, the non-parallel 

part of the code, is about 10s for the Yamabe equation and 5s for the f/fl steps on our SGI 

Origin 2000 with MIPS R10000 processors. 
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FIG. 5: Plots of the numerical violation of the constraint (II/14d,a = x), which is the 
constraint with the largest violation, on a 50 3 time evolution grid for the data calculated 
from a 17 x 16 x 16 (dotted line), a 25 x 24 x 24 (dashed line, almost invisible, since very 
close to the abscissa), and a 33 x 32 x 32 (solid line) spectral grid. 



For the asymptotically Minkowski case we have chosen 



o = i 

2 
/ 



h 



ab 



k 



2 i 2 i 2 

x + y + z 



1 + fi 2 cos y 



\VL 2 ((cosx) 2 — (siny) s 



V 



±0 2 (cosz) 2 



\Vt 2 ((cos a;) 2 — (siny) 5 
l + fi 2 (sin^) 2 
\Vt 2 cos a; 

5 



COS Z. 



\VL 2 [cosz) 



4 ^ 2 ( 

|fi 2 cosx 
1 + iO 2 



(36a) 

(36b) 
(36c) 



Without the fl 2 terms, the choice of hab leads to data for Minkowski space. 
Figure |6| shows the average violation (11/30) of the constraints. In this example the viola- 
tion of the constraints is not dominated by a single constraint as in the asymptotically A3 
example. The class of constraints with the largest violations are (II/14j), then followed by 
(II/14d) and (II/14g). 

As in the asymptotically A3 example the violation of the constraints is rapidly decreasing 
with the number of gridpoints in the spectral grid. For figure |6| we used a 34 x 32 x 32 
spectral grid. The violation of the constraints drops by a factor of ~ 8, when we refine the 
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FIG. 6: Average violation of the con- 
straints on a 50 3 (dashed line) and a 100 3 
(solid line) time evolution grid for data 
calculated from a 34 x 32 x 32 spectral 
grid. 



time evolution grid, on which we numerically evaluate the constraints, from 50 3 to 100 3 . If 
the data were an exact solution of the constraints, this factor would be ~ 16. If the data 
were a bad solution of the constraints, this factor would be ~ 1. Therefore we conclude that 
the data must be very close to an exact solution of the constraints. 



IV. THE CASE OF MULTIPLE BLACK HOLE SPACETIMES 

In the previous two subsections we have described in detail how initial data for the conformal 
field equations can be constructed for the asymptotically Minkowski and the asymptotically 
A3 case. The cases described are spacetimes with gravitational wave content. Depending on 
the "size" of the data, the gravitational waves may disperse to null infinity or may interact 
to form black holes. For our understanding of the Einstein equation it is certainly very 
interesting to study those cases, but as models for sources of gravitational waves they are 
probably not the most interesting ones. Therefore, the method for calculating initial data 
as described would only be of limited use, if it could not be extended to the case of one or 
more black holes (see figures 2 and 3 of I). In the following we describe how one can build 
a code to calculate data describing multiple black holes. 

The case of one black hole with radiation (asymptotically Schwarzschild) is a straightfor- 
ward extension of the asymptotically Minkowski case, since we can also use polar coordi- 
nates. Without loss of generality we can choose Q such that the inner Si coincides with 
the coordinate sphere at > and the outer S Q coincides with the coordinate sphere at 
r Q > ?y Instead of using [—1, 1] x [0, tt[x [0, tt[ to cover the physical part of the grid, we use 
[rj,r ] x [0, 7r[x[0, 2ir[. Regularisation is even easier, since r = is not part of E to , there is 
only an axis singularity but not an origin singularity to deal with. 

When extending to the time evolution grid inside Si more care is needed than in the asymp- 
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totically Minkowski case, since we extend onto the r = coordinate singularity of the time 
evolution grid. When extending towards r = the limit must not depend on the direction. 

The case of two black holes is significantly more complex. Since we do not see how to cover 
S( with a coordinate system for which we can also provide a spectral basis and fast spectral 
transformation algorithms, we suggest to use domain decomposition techniques. 
In figure [7| we give our three coordinate patches. One boundary of each patch coincides with 




FIG. 7: Patches for pseudo-spectral grids to cover the initial slice of a 
two black hole spacetime. Dashed lines mark the interior boundaries 
on which boundary values are read off from the other patches. 



a part of S, the other, the dashed lines, in future called Q, lies inside E to and inside at least 
one of the other patches. If we knew boundary values on Q, we would have three copies of 
the asymptotically Schwarzschild case — the elliptical shape of patch 1 is easily dealt with 
by a straightforward coordinate transformation. In the case of the f /Q divider we do know 
the boundary values at Q, namely f/fl, for the Yamabe equation we do not the boundary 
values at Q. 

To provide boundary values on Q we can use the Schwarz alternating procedure as described 
in section 6.4.1 of ||16|| : When solving the equation in one patch we read off the boundary 



values on Q from the grids covering the other patches. By iteratively solving the equation on 
all patches we calculate a solution on E to . The book ]16| contains a proof, that the Schwarz 
alternating procedure converges for the Laplace equation. The convergence rate depends on 
the overlap of the patches, the larger the overlap the better. In the overlap regions of our 
patches equation (g) is a normal elliptic equation, the principal part does not degenerate 
there. Therefore, it should in principle be possible to also prove convergence of the Schwarz 
alternating procedure for equation (Q). 
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V. CONCLUSION 

In this paper we have described a highly efficient and accurate scheme to calculate data for 
the conformal field equations without making any symmetry assumptions. The data are 
specified by giving a boundary defining function fl and the six components of the conformal 
3- metric h a b. 

Acknowledgement 

I would like to thank H. Friedrich for his help and support and J. Frauendiener for many 
helpful discussions on pseudo-spectral methods. 

B. Schmidt, M. Weaver, and S. Husa should also be mentioned for many discussions. 
I also acknowledge K. Stiiben from the Gesellschaft fiir Mathematik und Datenverarbeitung, 
who put the Algebraic Multigrid Library AMG at my disposal, and M. Frigo and S. G. John- 
son who wrote FFTW and made it publically available for the scientific community. 



[1] P. Hiibner, Class. Quantum Grav. 16, 2145 (1999). 
[2] P. Hiibner, Class. Quantum Grav. 16, 2823 (1999). 

[3] P. Hiibner, Numerische und analytische Untersuchungen von (singuldren, ) asymptotisch 
flachen Raumzeiten mit konformen Techniken, Ph.D. thesis, Ludwig-Maximilians-Universitat 
Miinchen (1993). 

[4] P. Hiibner, Phys. Rev. D 53(2), 701 (1996). 

[5] J. Frauendiener, Phys. Rev. D 58(6), 064002/1 (1998). 

[6] J. Frauendiener, Phys. Rev. D 58(6), 064003/1 (1998). 

[7] J. Frauendiener, J. Comp. Appl. Math. 109, 475 (1999). 

[8] J. Frauendiener, Liv. Rev. in Relativity (2000). 

[9] L. Andersson, P. T. Chrusci'el, and H. Friedrich, Comm. Math. Phys. 149, 587 (1992). 
[10] H. Friedrich, Commun. Math. Phys. 91, 445 (1983). 
[11] P. Hiibner, Class. Quant. Grav. 12, 791 (1995). 
[12] P. Hiibner and M. Weaver, in preparation. 
[13] K. Stiiben, GMD report (1999). 

[14] J. W. Ruge and K. Stiiben, in Multigrid Methods, edited by S. F. McCormick, SIAM (Society 

for Industrial and Applied Mathematics, Philadelphia, 1987), pp. 73-130. 
[15] D. Braess, Finite Elemente (Springer, Berlin, 1997). 

[16] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, vol. 23 
of Springer Series in Computational Mathematics (Springer, Berlin, 1997), 2nd ed. 

[17] S. Bonazzola, J. Frieben, E. Gourgoulhon, and J. A. Marck, gr-qc 9604029, 1 (1996). 

[18] M. Frigo and S. G. Johnson, in 1998 ICASSP conference proceedings, Vol. 3 (ICASSP, 1998), 
p. 1381. 

[19] C. Van Loan, Computational Frameworks for the Fast Fourier Transform, vol. 10 of Frontiers 

in Applied Mathematics (SIAM, Philadelphia, 1992). 
[20] S. Bonazzola and J.-A. Marck, J. Comp. Phys. 87, 201 (1990). 
[21] L. Andersson and P. T. Chrusciel, Dissertationes Mathematice pp. 1-100 (1996). 



21 



The choice of the sign is arbritrary, since the conformal and the physical metric are connected 
by the square of the conformal factor. 

We describe the 3D code here. We also have a 2D code, where we made use of the simplifications 
implied by an hypersurface orthogonal Killing vector d z a . From the description of the 3D code, 
the 2D code should be obvious. 



